## Author: Carl Lieberman
## Date: 2022-05-02

##This script produces individual output files that can be combined by a separate program.

##devtools::install_github("kolesarm/RDHonest@6de8ae146cd88280c86dbeaa7b136506da037589")

##Load required libraries:
library("Hmisc")
library("RDHonest")
library("haven")
library("foreign")
library("dplyr")

## Initialize directories
main <-"\\\\rschfs1x\\userRS\\a-e\\clieberman_RS\\Documents\\rd graph\\replication_final_public_test5"
dat <- "\\\\rschfs1x\\userRS\\a-e\\clieberman_RS\\Documents\\rd graph\\replication_final_public_test5\\data"
sim_dat <- paste0(main, "\\simulated_data")

## Create empty output file to populate
estimates <- data.frame()

## True derivatives
deriv <- read_dta(paste0(dat, "\\rdd_dgp_true_second_derivatives.dta"))

for (paper in c("ADGP1", "BDGP2", "CDGP3", "DDGP4", "EDGP5", "FDGP6", "GDGP7", "HDGP8", "IDGP9", "JDGP10", "KDGP11")){
    print(paste("The DGP is", paper))
    for (jump in c(-15, -9, -54, -324, -1944, 0, 1944, 324, 54, 9, 15 )){
        print(paste("This discontinuity is", jump))
        for (seed in c(1, 2, 3, 4, 5, 6, 7, 8)){
            print(paste("This seed is", seed))

            ## Load paper microdata
            dgpdata <- read_dta(paste0(sim_dat, "\\", paper, "_", jump, "_", seed, ".dta"))
            dgpdatarot <- RDData(data.frame(dgpdata$y, dgpdata$x), cutoff = 0)
            
            ## Bounds
            bound_taylor_true <- deriv$deriv_taylor[deriv$Paper == paper]
            bound_rot <- NPR_MROT.fit(dgpdatarot)

            ##Microdata results:
            rdh_taylor_true <- RDHonest(y ~ x, data = dgpdata, kern = "uniform",
                                        M = bound_taylor_true , opt.criterion = "MSE",
                                        sclass = "T", alpha=0.05)
            rdh_taylor_true_reject <- !(rdh_taylor_true$estimate - rdh_taylor_true$hl < 0 &
                                        rdh_taylor_true$estimate + rdh_taylor_true$hl > 0)
            rdh_taylor_rot <- RDHonest(y ~ x, data = dgpdata, kern = "uniform",
                                          M = bound_rot , opt.criterion = "MSE",
                                          sclass = "T", alpha=0.05)
            rdh_taylor_rot_reject <-
                !(rdh_taylor_rot$estimate - rdh_taylor_rot$hl < 0 &
                      rdh_taylor_rot$estimate + rdh_taylor_rot$hl > 0)
            
            results <- data.frame(paste0(paper, "_", jump, "_", seed),
                                  jump,
                                  rdh_taylor_true$estimate,
                                  rdh_taylor_true$sd,
                                  rdh_taylor_true$estimate - rdh_taylor_true$hl,
                                  rdh_taylor_true$estimate + rdh_taylor_true$hl,
                                  rdh_taylor_true_reject,
                                  bound_taylor_true,
                                  rdh_taylor_rot$estimate,
                                  rdh_taylor_rot$sd,
                                  rdh_taylor_rot$estimate - rdh_taylor_rot$hl,
                                  rdh_taylor_rot$estimate + rdh_taylor_rot$hl,
                                  rdh_taylor_rot_reject,
                                  bound_rot)
            
            colnames(results) <- c("file",
                                   "disc",
                                   "rdh_taylor_true_estimate",
                                   "rdh_taylor_true_sd",
                                   "rdh_taylor_true_lower",
                                   "rdh_taylor_true_upper",
                                   "rdh_taylor_true_reject",
                                   "rdh_taylor_true_bound",
                                   "rdh_taylor_rot_estimate",
                                   "rdh_taylor_rot_sd",
                                   "rdh_taylor_rot_lower",
                                   "rdh_taylor_rot_upper",
                                   "rdh_taylor_rot_reject",
                                   "rdh_taylor_rot_bound")

            estimates <- rbind(estimates, results)
            
            }
        }
    }

## Make sure paper name is a string, not a factor
estimates$file <- as.character(estimates$file)

var_labels = c(file="Dataset",
               disc="Discontinuity magnitude indicator",
               rdh_taylor_true_estimate="RDHonest discontinuity point estimate (w/ true deriv bound)",
               rdh_taylor_true_sd="RDHonest SD of discontinuity point estimate (w/ true deriv bound)",
               rdh_taylor_true_lower="RDHonest 95% CI lower bound (w/ true deriv bound)",
               rdh_taylor_true_upper="RDHonest 95% CI upper bound (w/ true deriv bound)",
               rdh_taylor_true_reject="RDHonest reject null indicator (p=0.05, w/ true deriv bound)",
               rdh_taylor_true_bound="True Taylor bound on DGP's 2nd deriv (requires omniscience)",
               rdh_taylor_rot_estimate="RDHonest discontinuity point estimate (w/ rule-of-thumb deriv bound)",
               rdh_taylor_rot_sd="RDHonest SD of discontinuity point estimate (w/ ROT deriv bound)",
               rdh_taylor_rot_lower="RDHonest 95% CI lower bound (w/ ROT deriv bound)",
               rdh_taylor_rot_upper="RDHonest 95% CI upper bound (w/ ROT deriv bound)",
               rdh_taylor_rot_reject="RDHonest reject null indicator (p=0.05, w/ ROT deriv bound",
               rdh_taylor_rot_bound="Rule-of-thumb Taylor bound on DGP's 2nd deriv")
label(estimates) = as.list(var_labels[match(names(estimates), names(var_labels))])

write.dta(estimates, paste0(dat, "\\rdhonest_estimates.dta"), version = 10)
